7393 Exam 2 R Code Andrei Matveev

library(data.table)
library(truncnorm)
library(coda)
devtools::install_github("stan-dev/loo")
library(loo)
library(rstan) 
library(shinystan)
library(invgamma)
library(bayesplot)
library(posterior)

devtools::install_github("jfrench/bayesutils")

setwd('/Users/AM/Documents/_CU Masters/2020 fall Bayesian_7393/code/Bayesian_Statistics_Class_Code/Exam_code')
#df1 = load("diamonds_simple.rda")
#rm(df1, diamonds_simple)
df = bayesutils::diamonds_simple

df$lprice = log(df$price)
df$lcarat = log(df$carat)

df_covid = bayesutils::covid_dec4

lq_theta_y = function(sigmaSQ, beta0, beta1, lpr = df$lprice, lcar = df$lcarat) {
  ld = dnorm(beta0, 0, 100, log = TRUE) + 
    dnorm(beta1, 0, 100, log = TRUE) + 
    dinvgamma(sigmaSQ, shape = 0.01, rate = 0.01, log = TRUE) +
    sum(dnorm(lpr, mean = (beta0 + beta1 * lcar), sd = sqrt(sigmaSQ), log = TRUE))
  return(ld)
}

mh = function(B, theta_start) {
  theta = array(0, c((B+1), 3), dimnames = list(c(), c("sigmaSQ", "beta0", "beta1")))
  
  theta[1,1] = theta_start[1]
  theta[1,2] = theta_start[2]
  theta[1,3] = theta_start[3]

  for (i in 2:dim(theta)[1]) {
    
    ###  step for sigmaSQ
    beta0_star = theta[(i-1),2]
    beta1_star = theta[(i-1),3]
    
    sigmaSQ_star = rtruncnorm(n = 1, a=0, b=Inf, mean = theta[(i-1),1], sd = 0.1)

    num_logr = lq_theta_y(sigmaSQ = sigmaSQ_star, beta0 = beta0_star, beta1 = beta1_star) -
      log(dtruncnorm(x = sigmaSQ_star, a=0, b=Inf, mean = theta[(i-1),1], sd = 0.1))
    den_logr = lq_theta_y(sigmaSQ = theta[i-1, 1], beta0 = beta0_star, beta1 = beta1_star) -
      log(dtruncnorm(x = theta[i-1, 1], a=0, b=Inf, mean = sigmaSQ_star, sd = 0.1)) 
## !!! is the mean correct? Should  we use in den_logr mean = sigmaSQ_star or mean = theta[i-1, 1]? 
# check c-componentwise-mh v1 line 68, 72
    logr = num_logr - den_logr
    if (log(runif(1)) <= min(logr, 0)) {
      theta[i,1] = sigmaSQ_star
    } else {
      theta[i,1] = theta[(i - 1), 1]
    }
    
    ###  step for beta0
    beta1_star = theta[(i-1),3] # it is the repeated code, but I need it to keep the interpretability
    sigmaSQ_star = theta[i,1] # update sigmaSQ_star after the Gibbs step for sigmaSQ
    
    beta0_star = rnorm(1, theta[(i-1),2], 0.1) # !!! check the parametrization (0.1 or 0.1^2)
    
    num_logr = lq_theta_y(sigmaSQ = sigmaSQ_star, beta0 = beta0_star, beta1 = beta1_star) -
      dnorm(x = beta0_star, mean = theta[(i-1),2], sd = 0.1, log = TRUE)
    den_logr = lq_theta_y(sigmaSQ = sigmaSQ_star, beta0 = theta[i-1, 2], beta1 = beta1_star) -
      dnorm(x = theta[(i-1),2], mean = beta0_star, sd = 0.1, log = TRUE)
## !!! is the mean correct? Should  we use in den_logr mean = sigmaSQ_star or mean = theta[i-1, 1]? 
# check c-componentwise-mh v1 line 68, 72
    logr = num_logr - den_logr
    if (log(runif(1)) <= min(logr, 0)) {
      theta[i,2] = beta0_star
    } else {
      theta[i,2] = theta[(i - 1), 2]
    }
    
    ###  step for beta1
    sigmaSQ_star = theta[i,1] # it is the repeated code, but I need it to keep the interpretability
    beta0_star = theta[i,2] # update beta0_star after the Gibbs step for beta0

    beta1_star = rnorm(1, theta[(i-1),3], 0.1) # !!! check the parametrization (0.1 or 0.1^2)
    
    num_logr = lq_theta_y(sigmaSQ = sigmaSQ_star, beta0 = beta0_star, beta1 = beta1_star) -
      dnorm(x = beta1_star, mean = theta[(i-1),3], sd = 0.1, log = TRUE)
    den_logr = lq_theta_y(sigmaSQ = sigmaSQ_star, beta0 = beta0_star, beta1 = theta[i-1, 3]) -
      dnorm(x = theta[(i-1),3], mean = beta1_star, sd = 0.1, log = TRUE)
## !!! is the mean correct? Should  we use in den_logr mean = sigmaSQ_star or mean = theta[i-1, 1]? 
# check c-componentwise-mh v1 line 68, 72
    logr = num_logr - den_logr
    if (log(runif(1)) <= min(logr, 0)) {
      theta[i,3] = beta1_star
    } else {
      theta[i,3] = theta[(i - 1), 3]
    }
    
  }
  return(theta)
}

B = 10^5 
keep = (B/2 + 1):(B + 1)
chain1 = mh(B, theta_start = c(0.1, -1, -1))
chain2 = mh(B, theta_start = c(0.3, 0, 0))
chain3 = mh(B, theta_start = c(0.5, -1, 1))
chain4 = mh(B, theta_start = c(0.2, 1, 1))

mc = mcmc.list(mcmc(chain1[keep,]), mcmc(chain2[keep,]),
                 mcmc(chain3[keep,]), mcmc(chain4[keep,]))
summary(mc)

Iterations = 1:50001
Thinning interval = 1 
Number of chains = 4 
Sample size per chain = 50001 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

          Mean      SD  Naive SE Time-series SE
sigmaSQ 0.1152 0.00925 2.068e-05      6.977e-05
beta0   8.3798 0.02240 5.008e-05      1.647e-04
beta1   1.5064 0.03135 7.011e-05      2.141e-04

2. Quantiles for each variable:

           2.5%    25%    50%    75%  97.5%
sigmaSQ 0.09853 0.1088 0.1145 0.1211 0.1348
beta0   8.33546 8.3649 8.3797 8.3948 8.4234
beta1   1.44522 1.4854 1.5064 1.5276 1.5680

gelman.diag(mc, autoburnin = FALSE)

Potential scale reduction factors:

        Point est. Upper C.I.
sigmaSQ       1.03       1.08
beta0         1.00       1.01
beta1         1.01       1.03

Multivariate psrf

1.03

# extracting the MCMC samples from the soda_fit object
samples = extract(r_fit_A)
ncycles = length(samples[[1]])
T <- length(df$price)

# each row of yrep is a sample from the pp distribution
yrep = matrix(0, ncol = T, nrow = ncycles)
for (i in seq_len(T)) {
  mui = samples$beta0 + samples$beta1 * df$lcarat[i]
  yrep[, i] = rnorm(ncycles, mean = mui, sd = sqrt(samples$sigmaSQ))
}

# approximate posterior predictive density for an observation
# with the same covariates as observation 1
#plot(density(yrep[,1]))
color_scheme_set("red")

# posterior predictive check
y = df$lprice
ppc_hist(y, yrep[sample(1:ncycles, 8),]) + ggtitle("Model A: histogram comparing the response to 8 replicated data sets")

# scatterplot of y vs average yrep
ppc_scatter_avg(y, yrep) + ggtitle("Model A: scatterplot comparing the response to the average replicated data set")

ppc_dens_overlay(y, yrep = yrep[sample(1:ncycles, 100),]) + ggtitle("Model A: density plot comparing the density of the response to the densities of 100 replicated data sets")

ppc_error_scatter_avg(y, yrep = yrep) + ggtitle("Model A: scatter plot comparing the response to the average predictive error")

ppc_intervals(y, yrep = yrep)

ppc_error_scatter_avg_vs_x(y = y, yrep = yrep, x = df$lcarat)

ppc_intervals(y = y, yrep = yrep, x = df$lcarat) + ggtitle("Model A: Posterior predictive intervals to the observed data values") + xlab("log carat") + ylab("log price")

# extracting the MCMC samples from the soda_fit object
samples = extract(r_fit_B)
ncycles = length(samples[[1]])
T <- length(df$price)

# each row of yrep is a sample from the pp distribution
yrep = matrix(0, ncol = T, nrow = ncycles)
for (i in seq_len(T)) {
  mui = samples$beta0 + samples$beta1 * df$carat[i] + samples$beta2 * df$carat[i] * df$carat[i]
  yrep[, i] = rnorm(ncycles, mean = mui, sd = sqrt(samples$sigmaSQ))
}

# approximate posterior predictive density for an observation
# with the same covariates as observation 1
#plot(density(yrep[,1]))
color_scheme_set("green")

# posterior predictive check
y = df$price
ppc_hist(y, yrep[sample(1:ncycles, 8),]) + ggtitle("Model B: histogram comparing the response to 8 replicated data sets")

# scatterplot of y vs average yrep
ppc_scatter_avg(y, yrep) + ggtitle("Model B: scatterplot comparing the response to the average replicated data set")

ppc_dens_overlay(y, yrep = yrep[sample(1:ncycles, 100),]) + ggtitle("Model B: density plot comparing the density of the response to the densities of 100 replicated data sets")

ppc_error_scatter_avg(y, yrep = yrep) + ggtitle("Model B: scatter plot comparing the response to the average predictive error")

ppc_intervals(y = y, yrep = yrep, x = df$carat) + ggtitle("Model B: Posterior predictive intervals to the observed data values") + xlab("carat") + ylab("price")

samples = extract(r_fit_D)
ncycles = length(samples[[1]])
T <- length(df$price)
df$notic = as.integer(df$inclusions == "noticeable")

# each row of yrep is a sample from the pp distribution
yrep = matrix(0, ncol = T, nrow = ncycles)
for (i in seq_len(T)) {
  mui = samples$beta0 + samples$beta1 * df$lcarat[i] + samples$beta2 * df$notic[i]
  yrep[, i] = rlnorm(ncycles, mean = mui, sd = sqrt(samples$sigmaSQ))
}

color_scheme_set("blue")
ppc_intervals_grouped(y = df$price, yrep = yrep, x = df$carat, prob = 0.5, prob_outer = .95, group = df$notic) + ggtitle("Model D: Posterior predictive intervals to the observed data values") + xlab("lcarat") + ylab("price")


find_ep = function(yrep, min = 0.1, max = 2.4, step=0.1) {
  ep = array(0, dim=c((max-min+step)/step, 5), dimnames = list(c(), c("carat_min", "carat_max", "2.5%","Mean","97.5%")))
  n=0
  for (i in seq(min+step, max, by=step)) {
    n=n+1
    subset = ((df$carat>(i-step)) & (df$carat<=i))
    yrep_subset = yrep[,subset]
    ep[n,1] = i-step
    ep[n,2] = i
    q = quantile(yrep_subset, probs = c(0.025, 0.975))
    ep[n,3] = q[[1]]
    ep[n,4] = mean(yrep_subset)
    ep[n,5] = q[[2]]
  }
return(ep)
}

find_ep(yrep)
temp = ppc_intervals_data(y = df$price, yrep = yrep, x = df$lcarat, prob = 0.5, prob_outer = .95, group = df$notic)

```